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Abstract 

We present a lattice investigation of the heavy-Ught meson de- 
cay constants using Wilson light quarks and NRQCD heavy quarks, 
partially including the effects of dynamical sea quarks. We calcu- 
late the pseudoscalar and vector decay constants over a wide range in 
heavy quark mass and are able to perform a detailed analysis of heavy 
quark symmetry. We find consistency between the extrapolation of the 
NRQCD results and the static case, as expected. We find the slope of 
the decay constants with 1/M is significantly larger than naive expec- 
tations and the results of previous lattice calculations. For the first 
time we extract the non-perturbative coefficients of the slope arising 
from the 0{1/M) heavy quark interactions separately and show the 
kinetic energy of the heavy quark is dominant and responsible for the 
large slope. In addition, we find that significant systematic errors re- 
main in the decay constant extracted around the B meson mass due 
to truncating the NRQCD series at 0(1/M). We estimate the higher 
order contributions to fs are approximately 20%; roughly the same 
size as the systematic errors introduced by using the Wilson action 
for light quarks. 
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1 Introduction 



The decay constant of the B meson, Jb, is of fundamental importance to 
tests of the Standard Model because it is needed in the determination of the 
CKM matrix. While one of the simplest weak matrix elements to study, it 
is not very well known and dominates the uncertainty in extracting Vtd, one 
of the smallest CKM elements, from B — B mixing. The best chance of an 
experimental measurement of lies in the decay B^ ^ r+ + z/. However, 
this decay has yet to be observed, and appears in combination with Vuh, 
which is only known to within a factor of 2. 

Theoretical calculations using for example QCD sum rules within the 
framework of Heavy Quark Effective Theory (HQET) [Q, || and lattice field 
theory have been actively pursued over the last few years with the current 
world average standing at fs = 200 ± 20% MeV 0, ^. As a first principles 
approach, lattice simulations offer ultimately the most reliable calculation of 
fs- The initial problem of simulating a h quark with Compton wavelength of 
the order of the typical lattice spacing presently achievable was first solved 
by using the static approximation , where the heavy quark is treated as 

a static colour source within the heavy-light meson. Simulations using con- 
ventional relativistic methods around the D meson showed the corrections to 
the static approximation in this region, and hence also at the B meson, are 
much larger than naively expected. This has contributed to the sizable un- 
certainty remaining in /s, shown above. The determination of fs directly at 
Mb on the lattice is required, and this is now possible with the development 
of Non Relativistic QCD (NRQCD) 0, |^. Using this approach the error 
in the decay constant can be significantly reduced with the computational 
power currently available. Initial calculations have already been performed 



by Hashimoto [p| and Davies |10 



This analysis forms the second stage of a project applying NRQCD to 
heavy-light systems and studying heavy quark symmetry at the B meson. 
Our previous papers have dealt with the B meson spectrum in both the 
quenched approximation [|lT| and partially including the effects of sea quarks 
[p!^ . We now present our results for the decay constants of the heavy-light 
mesons from dynamical configurations. Note that the combination /^-Bs, 
where Bb is the mixing parameter for the B system, is the experimentally 
relevant quantity, and thus a calculation of Bb is also needed. 

The paper is organised as follows. In section || we describe the simulation 
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details including the systematic errors associated with the calculation. The 
meson operators required to compute the decay constant to 0(1/M) are set 
out in section |^ and we discuss computing the decay constant consistently 
to this order in terms of the evolution equation used. Next, we examine the 
predictions of heavy quark symmetry for the behaviour of the decay constant 
in the heavy quark mass limit. The contributions from the interactions of the 
heavy quark at 0(1/M) are analysed and combinations of decay constants 
which isolate these individual contributions are discussed. Section |^ illus- 
trates our results and fitting analysis with a subset of the data. The heavy 
quark mass dependence of the decay constants are investigated in section 
the individual non-perturbative 0(1/M) coefficients of the slope of the decay 
constant are extracted. Finally we present our conclusions in section ^. 



2 Simulation Details 

The simulations were performed using 100 16^ x 32 gauge configurations at 
j3 = 5.6 with two fiavours of staggered dynamical sea quarks with a bare 
quark mass of arrisea = 0.01. These configurations were generously made 



available by the HEMCGC collaboration; more details can be found in [|I3 
We fixed the configurations to the Coulomb gauge. The light quark prop- 
agators were generated using the Wilson fermion action, without an 0(a) 
improvement term, at two values of the hopping parameter, k = 0.1585 and 
0.1600. The former corresponds to a quark mass close to strange, where 
Us = 0.1577 from M^, while 0.1600 is somewhat lighter, with Kc = 0.1610. 

With only two values of light quark mass it is not possible to perform 
a trustworthy chiral extrapolation. In addition, the systematic errors asso- 
ciated with the light quarks, discussed below, affect the determination of 
the bare strange quark mass; Ks extracted from pseudoscalar mesons, using 
the lowest order chiral mass dependence to find the 'experimental mass' of 
the pure ss pseudoscalar, differs from Kg obtained from the vector meson (p. 
Thus, most of the results will be presented for ni = 0.1585. 

In this simulation we truncate the NRQCD series at 0{1/Mq), where Mq 
is the bare heavy quark mass, and the action takes the form: 

S = Q\Dt + Ho + 6H)Q (1) 
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where 

A(2) a ■ B 

Ho = and 6H = -cb . (2) 

° 2Mo 2Mo ^ ^ 

Tadpole improvement of the gauge hnks is used throughout i.e. they are 
divided by uq, where we use Uq = 0.867 measured from the plaquette, and 
the hyperfine coefficient is given the tree-level value cb = 1. We use the 
standard Clover-leaf operator for the B field. The heavy quark propagators 



were computed using the evolution equation [|14 



2n J "V 2n 



n 



Gt+1 = {^-<^mGt (t>l) (3) 

where n is the stabilising parameter. In the static limit this reduces to 

Gt+i - UlGt = 4,0- (4) 

We generated heavy quark propagators at 11 values of (aMo,n) corre- 
sponding to (0.8,5), (1.0,4), (1.2,3), (1.7,2), (2.0,2), (2.5,2), (3.0,2), (3.5,1), 
(4.0,1), (7.0,1) and (10.0,1), and the static limit. This roughly corresponds 
to a range of meson masses from Mb/3 to AMb and is sufficient for a reason- 
able investigation of heavy quark symmetry. However, it is not possible to 
simulate the D meson on this lattice using NRQCD; a larger lattice spacing 
is required, /3'"^=2<5.5 or /5"^=o<5.85. 

Details of the construction of meson operators, the smearing functions 



used and our fitting analysis can be found in [|I2| . As discussed in this refer- 
ence there are large systematic errors associated with this simulation arising 
from the use of Wilson light fermions. The emphasis of our analysis will be 
on determining the heavy quark mass dependence of the decay constants and 
extracting the nonperturbative coefficients of the heavy quark expansion to 
0(1/M). The systematic error arising from the light quark in these 0(Aq(^^) 
coefficients is naively estimated to be a relative error of 0{Kqcd<i) ~ 20% 
for this ensemble (where we take ah.QCD = ci-^v = 0.185 for these configura- 
tions). The truncation of the NRQCD series and the use of the tree- level value 
for the hyperfine coefficient introduce absolute errors of 0{Aqcd{-^qcd/ M^) 
and O^asAqcDi-^QCD/M)) respectively. Both correspond to approximately 
1% errors in the value and 10% errors in the slope of the decay constant at 
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Mb- The coefficients of the quadratic term in the heavy quark expansion 
cannot be correctly determined in this analysis. The size of the systematic 
errors induced by the light quark action are indicated by the large uncertainty 
in the lattice spacing; we use = 1.8—2.4 GeV. This gives oMq = 2.4—1.7 
as corresponding to the bare b quark mass. We are currently repeating our 
analysis using tadpole- improved Clover light fermions ||r5| . 



3 Meson Operators 

The pseudoscalar decay constant is defined by 

{0\Ao\PS) = fpsMps, (5) 

where Aq = q'j^'joh at tree level and q and h represent 4-component light and 
heavy quark fields respectively. In this convention the experimental value of 
/tt is 131 MeV. Similarly, the vector decay constant is given by 

{0\V,.m = e,fvMv,' (6) 

where Vi = q'jih at tree level. The calculation of these matrix elements, and 
thus decay constants, to 0(1/M) in NRQCD requires several operators. At 
tree level these can be obtained by relating h of full QCD to the 2-component 
NRQCD fields, Q and Q for heavy quarks and heavy antiquarks respectively, 
using the Foldy-Wouthuysen transformation. To 0(1/M), 

.(1 .-^-il (7, 

With the convention that our meson states be built out of a heavy quark 
and a light antiquark, the Q field will not contribute to the above matrix 
elements through 0(1/M), and we obtain 

(0| qT^'^h \P)qc^ = - (0| qlT^'^Q \P)^nQCD+{0\ qlT^'^ {j]^^ ^ \P)nrqcd 

(8) 




■^This is a more convenient definition than the more commonly used {0\qjih\Vi) 
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where r*^'^-' = 7570 and T^"^^ = 1 for the axial current, and T^^^ = "fi and 
r^^) = iai for the vector current. qi2 and qu denote the upper and lower two 
components of q respectively. Our Euclidean space 7 matrices are chosen as, 



70= : " , , „ . (9) 





Beyond tree level all operators with the same quantum numbers can mix 
under renormalisation. Thus at 0{a/M) the basis of operators for the axial 
current, is given by f\ 

o, = -,lQ o, = 4.'^Q f3 = ^i^^0 (10) 

Using translation invariance, the zero momentum currents involving O2 and 
O3 are identical on the lattice and comparing with equation ^ it is sufficient 
to calculate matrix elements of the tree level operators Oi and 02- Similarly, 
the basis of operators for the vector current, 

Oi = -iqla.Q O, ^ q^a^^^^Q O, ^ ql,_^a,Q (11) 

is reduced, for the numerical computation, to Oi, O2 and C2— £^4, where a 
matrix properties have also been used. We computed the pseudoscalar and 
vector matrix elements corresponding to this minimal set of operators by 
inserting a smeared interpolating field with the same quantum numbers as 
Oi at the source and all combinations of local currents at the sink for each 
meson. 

The perturbative calculation of the renormalisation factors needed to re- 
late the matrix elements in NRQCD to the currents in full QCD has not yet 
been completed [|l6l. Thus, we use the static renormalisation factors and 



are restricted to the tree level 0(1/M) corrections; we omit O/^ for the vector 
meson. Since this correction is of 0{as^QCD{.^QCDl M)) it is roughly the 
same size as that induced by unknown perturbative matching corrections in 

"^In the limit of zero light quark mass. 
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the coefficient cb, which are not included in the analysis. The renormalisa- 
tion factors in NRQCD depend on the heavy quark mass and hence using the 
static value introduces an 0{as) ~ 20% error into the 0(1/ M) coefficients 
extracted from our results. 

It is possible that the power divergences which appear in NRQCD may 
lead to significant non-perturbative contributions to the renormalisation fac- 



tors 181. No clear indication of 'renormalon' effects has been found in other 



quantities calculated perturbatively in NRQCD |]T4|, however, we intend to 



compute the renormalisation factors non-perturbatively to check this. 

Strictly, our analysis is not fully consistent to 0(1/M). Note that the 
evolution equation given in equation ^ induces an 0{a{AQCD/M)) system- 
atic error in the currents. This is because it does not apply the relativistic 
correction operator (1 — 6H) on the initial timeslice. Since all subsequent 
timeslices include this operator the transfer matrix is consistent to 0{1/M). 
Thus, it is only in the amplitude of the meson correlator, i.e. the decay 
constant, that is missing this 0{a{AQCD/M)) contribution. A better but 
computationally more costly evolution equation is given by 

0,^. . (13) 

for all timeslices. Since our use of the Wilson action for light quarks also 
induces errors of this order of magnitude, we used equation |^ in this project. 
In a comparison of these evolution equations using Clover fermions for the 
light quarks ||l5l we found the insertion of (1 — 6H) on the source timeslice 



introduces ~ 4% correction to (0| Oi \PS) j^j^Qfj^ and ~ 3% correction to 
(0| O2 \PS) 

NRQCD around Mf,. Considering the systematic errors inherent in 
this calculation and our statistical accuracy of ~ 2%, the omission of this 
correction will not significantly affect our findings. 



4 Heavy Quark Symmetry 

In the heavy quark limit the decay constant can be parameterised in terms 
of an expansion in 1/M: 

/VM = (/VMr(l + | + o(^)). (14) 
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The coefficient cp is determined by nonperturbative contributions arising 
from the hyperfine interaction (Ghyp) and the kinetic energy of the heavy 
quark {Gkm), which appear in the NRQCD action, and the corrections to 
the current arising from the Foldy-Wouthuysen transformation. From first 
order perturbation theory in 1/M about the static hmit: 

Cp = Gkin + '^^MGhyp + dMGcorr/Q (15) 

where (Im = 3 and —1 for pseudoscalar and vector mesons respectively. The 
integer factors in this equation are introduced for convenience when compar- 
ing with HQET. 

(0|gtr(2)Q|p)^Gfc,„ = 

(0| / rf2/T{gtr(^)Q(0), QH^)Q{y)} l^)oo , (i6) 
{0\q^T^^^Q\P)^2dMGHyp = 

(0| / ciyT{gtr(^)g(o),gt( -^^^^ '^ )Q{y)} |P)^(i7) 



gtr(2)Q|P)^rf,,G,o,,/6 = (0|gtr(2)^Q|P) 



oo ' 



where \P)^ represents the meson in the hmit of infinite heavy quark mass, 
and equations |l6| - p!^ are tree level expressions. Note that the kinetic and 
hyperfine terms contribute to the decay constant through the correction to 
the meson wavef unction. 

Equation |15| suggests that the individual contribution from each 0{1/M) 
term can be obtained separately by taking appropriate combinations of the 
pseudoscalar and vector decay constants, {f\fM)ps and {fVM)v respec- 
tively. For clarity below, we denote 



{fVM)Ti- = -^jO\Or\PS)^j,QaD (19) 

1 



LPS 

lP5 



6{fVM)ps = ^j^moP\PS)^^Q^,^ (20) 



ifVM)fs = if^^)7s"'^'^ + S{fVM)ps (21) 
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expansion 


slope 


ifVM) 


(/VM)- (1 + c'p/M) 




( / ) fs'"''' / ( / v^) T""'''' 


1 + 4/M 


Cp = SGfiyp 


{fVM)%/{fVMy^' 


1 + c'^/M 


Cp = SGfiyp + 2GcQrr/3 



Table 1: Expansions about the static limit to 0(1/M) of combinations of 
the pseudoscalar and vector decay constant, with and without the current 
corrections. 



with similar quantities defined for the vector particle by replacing PS with 
V. We also denote the spin-average of the decay constants, 



(fVM) = ((/VM)^™ + 3(/VM)™") /4. (22) 

Table |I| displays our notations for the expansions in 1/M and the corre- 
sponding coefficients of (/a/M) and the ratio of the pseudoscalar to vector 
decay constants, with and without the current correction. It is easy to see 
from this that both Gkin and Ghyp can be extracted in a simple way from Cp 
and Cp respectively; Gcorr must be obtained from the combination Cp — Cp. 
Since we compute (0| O2 \P) separately, however, Gcorr can be more accu- 
rately found by examining the current correction in the infinite heavy quark 
mass limit: 



_ [2Mo(O|02|P)]oo _ J. 2Mo6{fVM)ps , . 

{0\O,\P)^ Mo™c» (/v/M)~ • ^ ^ 



HQET provides an analogous decomposition of these coefficients |I9| 
However, since HQET constructs an effective theory in terms of the heavy 
quark pole mass, as opposed to the bare quark mass in NRQCD, only physical 
combinations of G's, i.e. cp, Cp and Cp = Gkin, can be compared. Note that 
while these coefficients are expected to be of 0{Aqcd) ~ 200—500 MeV (and 
thus the corrections to the static limit at Mp 0{Aqcd/Mb) ~ 10%); the 
unphysical quantities Ghyp and Gcorr need not be of this order. 



In HQET Gcorr is related to the meson binding energy, Gcorr = —A [|1£ 
This follows from the observation 

{0\tdc.{qr('^h)\P)^ = {MM-MQ)v^{0\qr^'^h\P)^, (24) 
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= Av^{0\ qT('^h\P)^, (25) 

where Va is the velocity of the heavy quark, and using the equations of 
motion of the heavy and hght quark. Note that in NRQCD the R.H.S. 
of equation becomes E'^^Va (0| qV^^^h \P)^ and following the arguments 
of [19| the analogous expression in NRQCD can be obtained: 



Gcorr — ^sim- (26) 

At finite heavy quark mass this relation is broken by terms of 0{1/M). Sim- 
ilarly, on the lattice, discretisation modifies and the operators in equa- 
tion |23| must be renormalised. 



5 Results 

5.1 Fitting results for the decay constants 

To begin with we extract the decay constant without the correction to the 
current, i.e. from (0| Oi \P). We performed a single exponential vector fit to 
C{l,loi) and C(l,l) for the pseudoscalar and vector mesons. C(l,/oJ de- 
notes a correlation function constructed at the source from an interpolating 
field with appropriate quantum numbers for the vector or pseudoscalar state 
and the heavy quark smeared with a ground S-state smearing function (de- 
scribed in |T^); at the sink the local current Oi is inserted. In the limit of 



large values of t 

C(l,/oJ = Zi^^Z.e-''-', (27) 
C(l,l) = Z^e"^— *. (28) 

By taking a bootstrap ratio of the amplitudes of these correlators the decay 
constant can be obtained: 



V2Zi,^ = (0| \P)r,nQCD = (/v^)"- (29) 

To illustrate our results we use the data at aMo = 1.0 and ni = 0.1585. 
The effective masses of the pseudoscalar C{l,loi) and C(l,l) correlators, 
presented in figure |l|, show a clear signal out to approximately timeslice 
20. A one exponential {n^xp = 1) simultaneous fit to both correlators that 
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satisfies the fit criteria > 0.1 is possible from an initial timeslice (tmm) of 
4 out to timeslice 20 (tmax)- The variation of Esim with the fit range is ~ 2a 
for tmin^4: and this is taken to be the fitting error in our analysis. A final 
fitting range of 7 — 20 is chosen, for which Q is judged to be roughly stable 
to increases in tmin- 

To estimate the remaining excited state contribution to the fitting param- 
eters extracted the results are compared in table ^ with a multi-exponential 
analysis using a vector and matrix of smeared correlators [|12]. The agree- 
ment, to within the fitting error, of the energies and amplitudes provides 
confidence that we have isolated the ground state. Since the overlaps of both 
C{l,loi) and C(l, 1) with the first excited state are small, Uexp = 2 fits could 
not be performed successfully to these correlators. 



Data 


^exp 


fit range 


Q 


n F^^ 






11,11 


1 


7-20 


0.7 


0.487(2) 


1.09(2) 


10.1(1) 


11,21 


2 


3-20 


0.3 


0.484(3) 


1.02(3) 




11,12,21,22 


2 


3-20 


0.4 


0.485(2) 




10.0(1) 



Table 2: The fit parameters, ground state energies and amplitudes extracted 
from vector and matrix fits to the pseudoscalar meson correlators (without 
the current correction) for oMq = 1.0 and = 0.1585. 



In the static case, it is more difficult to obtain a satisfactory fit. From fig- 
ure 1^ the signal for the C(l, 1) dies away as C(l, /^J appears to be reaching 
a plateau. While a Uexp = 1 fit to both correlators is possible from tmin = 4 
for Q > 0.1, figure | shows the ground state energy is not stable to changes 
in the fitting range. In addition, it was not possible to bootstrap the fits for 
tmin < 9. Thus, we choose 9 — 11 for which Egim is consistent with the results 
from different fitting ranges. Agreement is found when comparing with the 
energy and amplitudes obtained from fits to other combinations of correla- 
tors, detailed in table ^. However, a better tuned ground state smearing 
function would enable the fit parameters to be extracted with greater con- 



fidence. As noted in reference [12|, even for the heaviest finite heavy quark 



^Defined as the one-standard deviation confidence level. 
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0.4 I < < < < ' < < < < ' < ^ 0.47 

10 20 5 10 15 

Tlmeslice t , 



1.0 




5 10 15 



Figure 1: The effective masses of the pseudoscalar C(l,ZoJ (circles) and 
C(l,l) (diamonds) correlators, and the corresponding ground state energy 
extracted as a function of tmin from a Uexp = 1 vector fit for gMq = 1.0 and 
Ki — 0.1585. The effective mass of the C(l, /oj correlator is offset by +0.2. 
Also shown are the values of Q for the vector fit; is only presented for 
those values of tmin which satisfy Q > 0.1. tmax is fixed at 20. 
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mass value, Mq = 10.0 shown in figure ^ NRQCD does not suffer from tlie 
same signal /noise problems. 



Data 


^exp 


fit range 


Q 








11,11 


1 


9-11 


0.65 


0.529(13) 


2.07(20) 


10.7(5) 


11,21 


1 


8-15 


0.2 


0.530(9) 


2.1(1) 




11 


1 


3-11 


0.6 


0.527(6) 




10.4(2) 



Table 3: The fit parameters, ground state energies and amplitudes extracted 
from vector fits to the pseudoscalar meson correlators for aMo = oo and 
Ki = 0.1585. 



The heavy quark is expected to be almost a spectator in the heavy-light 
meson and so the wavefunction changes very little with Mq or the spin content 
of the meson. We find it is sufficient to use the same fitting range, 7—20, for 
the pseudoscalar mesons for all finite values of the heavy quark mass. This 
fitting range is also used for the vector mesons, apart from Mq = 0.8—1.2 
where it was necessary to use 9—20. The corresponding values of Esim and 
the pseudoscalar and vector decay constants (without the current correction) 
are given in table ^ for ni = 0.1585. The pseudoscalar meson mass, Mps, 
also shown in the table, is calculated using the mass shifts, Mps—Egim, given 



in reference [0. Repeating the analysis for the results at ki = 0.1600, we 
found 7 — 20 is the optimal fitting range for all Mq. 

The current correction is extracted separately by taking the jackknife 
ratio of the pseudoscalar C{l,lo2) C{l,loi) correlators. In the limit of 
large times, the ratio tends to a constant: 



C(l,/aJ S{fVM)ps 



C(l,/oJ (/Vm)™^' 



(30) 



This ratio is shown in figure^ for aMo = 1-0 and ki = 0.1585. While a plateau 
seems to be reached around timeslice 5, Q does not stabilise until tmm = 12. 
Thus, we choose a fitting range of 12 — 20; this was found to be optimal 
for all Mq. In addition, picking an earlier tmm leads to difficulties when 
investigating the heavy quark mass dependence of the current correction. The 
statistical uncertainty becomes so small that terms above 0{1/M^) in the 



12 



10 ao 

Timeslice 



1.0 



0.6 



0.2 
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Figure 2: The effective masses of the pseudoscalar C(l,/oJ (circles) and 
C(l,l) (diamonds) correlators, and the corresponding ground state energy 
extracted as a function of tmin from a Uexp = 1 vector fit for aMg = oo and 
Ki = 0.1585. The effective mass of the C{l,loi) correlator is offset by +0.2. 
Also shown are the values of Q for the vector fit. tmax is fixed at 11. 
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0.6 



10 
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Figure 3: The effective masses of the pseudoscalar C{l,loJ (circles) and 
C(l, 1) (diamonds) meson correlators for uMq = 10.0 and ki = 0.1585. 

heavy quark expansion of this quantity are required, and it becomes difficult 
to extract the coefficients reliably. Using the previous results for {f\fM)ps 
we obtain the values of 5{f \fM)ps shown in table HJor k/ = 0.1585. It is 
clear that as 5{fy/M)ps is a ~ 15% correction to (/VM)p" around (~ 
2.0), it is essential to include the current correction when calculating fs 
and investigating its corresponding heavy quark mass dependence, especially 
when other sources of error from light quarks etc are under control. 

As stated in equations O and O, group theory relates the spin dependent 



contributions to the slopes of the vector and pseudoscalar decay constants in 
the continuum. In terms of the current correction: 



5(/v/M)p5/(/VM) 



lUncorr 
PS 



6{fVM)v/{fVM) 



1 uncorr 
V 



(31) 



On the lattice this relation only holds in the limit of infinite statistics. Fig- 



[C(l,lo^)/C{l,lo^)]ps 
[C{l,lo2)/C{l,lo-^)]v 



for aMo = 1.0 and ki = 0.1585. In the limit 



ure ^ shows 

of large times this ratio of correlators tends to the R.H.S. of equation ^ 
We find equation |31| is satisfied to within ~ l—2a and thus, we can further 
reduce the number of matrix elements that need to be calculated by using 



j^-^uncorr ^.j^^ vcctor curreut correction. 
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0.32 r 



0.280 r 



0.30 



0.28 



0.26 



0.274 



10 



20 



5 10 



15 



1.0 
0.8 
0.6 
0.4 
0.2 
0.0 

5 10 15 

tmln 

Figure 4: The jackknife ratio of the pseudoscalar C{l,lo2) and C{l,loi) 
correlators and the corresponding amphtude, 5{f\/M)ps/{f\fM)'ps°^^, as a 
function of tmin obtained from a constant fit for oMq = 1.0 and ki — 0.1585. 
The values of Q for the fits are also shown, tmax is fixed at 20. 
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aMo 


E ■ 




a3/2(/v^)«™ 




a3/2(/v^)™^ 


0.8 


0.475(3) 


1.36(2) 


0.180(3) 


-0.0605(13) 


0.207(5) 


1.0 


0.488(3) 


1.58(3) 


0.191(3) 


-0.0529(10) 


0.179(5) 


1.2 


0.497(3) 


1.76(2) 


0.199(3) 


-0.0470(9) 


0.188(5) 


1.7 


0.510(3) 


2.27(2) 


0.215(3) 


-0.0373(8) 


0.205(5) 


2.0 


0.514(3) 


2.58(2) 


0.224(3) 


-0.0337(8) 


0.216(5) 


2.5 


0.519(4) 


3.09(2) 


0.237(4) 


-0.0289(6) 


0.228(5) 


3.0 


0.522(3) 


3.52(3) 


0.247(5) 


-0.0255(8) 


0.239(5) 


3.5 


0.524(3) 


4.16(3) 


0.251(5) 


-0.0224(6) 


0.244(5) 


4.0 


0.525(3) 


4.64(4) 


0.258(5) 


-0.0204(6) 


0.252(5) 


7.0 


0.528(4) 


7.13(12) 


0.285(5) 


-0.0132(3) 


0.282(6) 


10.0 


0.527(4) 


9.95(12) 


0.300(6) 


-0.0099(3) 


0.297(8) 


oo 


0.528(5) 




0.341(37) 







Table 4: The pseudoscalar ground state energy, meson mass, decay constant 
and current correction, and the vector decay constant for all oMq and ki = 
0.1585. Note that the renormalisation factors are not included. 



We compare Ef^ with 2M, (0| O^^ \PS)^^qcd I (0| (^f ^ \PS)^nQCD ^ 
table |[ The 20% ~ 0{as) disagreement between these two quantities is 
approximately independent of heavy quark mass, and hence is most likely 
to be due to the omission of the lattice renormalisation factors rather than 
0(1/M) terms. Note that both quantities are roughly a factor of 2 larger 
than A i.e. quite different from the analogous quantities in HQET. 

The spin-average of the vector and pseudoscalar decay constants, as well 
as the ratio, with and without the current correction were computed. The 
results are given in table |^. In order to investigate the heavy quark mass 
dependence of these quantities, detailed in the next section, 100 bootstrap 
ensembles were generated for each correlated fit. The bootstrap averages 
of the fit parameters are consistent with the corresponding unbooted values 
given in the tables. In the static case we obtain, 

aE'^i^ = 0.559(37) a^ZiZi = 2.82(88) a^/^Zi = 12.4(1.7). (32) 

These values are consistent, but with much larger errors than those obtained 
from unbooted data. This underscores the signal /noise problems and the 
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Figure 5: The jackknife ratio of C(l, /o2)/C(l, /^J for the vector and pseu- 
doscalar mesons for oMq = 1.0 and ni = 0.1585. 



rough nature of our determinations of (/vM)'**"**'^. The corresponding esti- 
mate of {f\fMY*°-* = 0.4(1) is so uncertain as not to be useful Since the pur- 
pose of our static result is to compare with the more accurate extrapolation 
to the infinite mass limit of the NRQCD results, a rough value is sufficient 
and we use a value consistent with the unbooted fits, (/v^M)"^*"* = 0.341(37), 
which has smaller error bars. 

While the one loop correction to the decay constant is not included in this 
analysis, the size of the vector matrix element (0| O4 \ V) j^j^q(jj-, is shown in 
figure H as a percentage of (0| Oi |^)7Vii;QCD ^'^^ ~ ^-^ ~ 0.1585. 

This is roughly the same size as the tree level matrix element. 



5(/VM)v/(/VM) 



^uncorr 



M)ps/{f 



]\ r\uncorr 



0.28/3. (33) 



Since the one loop contributions should be suppressed by a factor of 0{as 
20%, the contribution to the vector decay constant is around 2% of 
(0| Oi \V)^j^Qfjj^. This is of order ~ 2a in the statistical errors of {f\fM) 
at this mass. The one loop correction to the pseudoscalar decay constant is 
similarly small. 



uncorr 
V 
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(lE • 




U.o 


()Aio[o ) 


U.t)o / (^zj 


1.0 


0.488(3) 


0.554(3) 


1.2 


0.497(3) 


0.567(3) 


1.7 


0.510(3) 


0.589(3) 


2.0 


0.514(3) 


0.600(3) 


2.5 


0.519(4) 


0.612(4) 


3.0 


0.522(3) 


0.621(4) 


3.5 


0.524(3) 


0.627(4) 


4.0 


0.525(3) 


0.631(5) 


7.0 


0.528(4) 


0.649(6) 


10.0 


0.527(4) 


0.657(8) 



Table 5: A comparison of the pseudoscalar ground S-state energy and 
{0\O2\PS)^^QCj^/ {0\Oi\PS)^j^QCjj for all oMq and ki = 0.1585. Note 
that the renormalisation factors are not included. 
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Figure 6: The jackknife ratio of C{1, 104,) /C{1, lo^) for the vector meson at 
aMo = 1.0 and ki = 0.1585. 
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aMo 




/ {fV^'^)pS \uncorr 


/ yj^j m )ps \tot 


0.8 


0.172(5) 


1.070(18) 


0.641(11) 


1.0 


0.182(5) 


1.065(19) 


0.706(13) 


1.2 


0.190(5) 


1.059(18) 


0.751(13) 


1.7 


0.208(3) 


1.047(7) 


0.819(6) 


2.0 


0.217(5) 


1.042(7) 


0.845(6) 


2.5 


0.229(5) 


1.037(7) 


0.875(5) 


3.0 


0.241(5) 


1.032(7) 


0.895(5) 


3.5 


0.246(5) 


1.028(6) 


0.910(5) 


4.0 


0.253(5) 


1.025(7) 


0.920(6) 


7.0 


0.282(6) 


1.013(5) 


0.952(4) 


10.0 


0.298(8) 


1.007(4) 


0.964(5) 



Table 6: The spin-averaged decay constant, and the ratio of the pseudoscalar 
to vector decay constant, both with and without the current correction for all 
aMo and k; = 0.1585. Note that the renormalisation factors are not included. 



5.2 Heavy quark mass dependence of the decay con- 
stants 

The decay constants have been computed over a wide range of heavy quark 
mass and this is ideal for investigating heavy quark symmetry. Figure |^ 
presents the results for (/v^M)p5 as a function of 1/Mps for ki = 0.1585. 
The static result is also included in the plot. A correlated fit to the data was 
performed using the functional form: 

/Vm = Co + ^ + ^ + ^. (34) 

The fitting range was varied keeping the starting point fixed at the heaviest 
data point, oMq = 10.0. Beginning with a linear function, the fitting pro- 
cedure was repeated adding quadratic and cubic terms. In general, in order 
to have confidence in the value for Ci a quadratic fit is required. Note that 
in reference |^2| we found when performing correlated fits to the results as a 
function of 1/M the values of n used in the evolution equation for oMq = 3.5 
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0.3 
"ffl 0.2 



0.1 

0.0 0.2 0.4 0.6 0.1 

l/(aMps) 

Figure 7: a^/^{f\fM)% vs l/{aMps) for ki = 0.1585. The solid line indicates 
a quadratic fit to the heaviest six points and l/laMs) ~ 0.4. The static 
data point, shown as a burst, is not included in the fit. Note that the 
renormalisation factors are not included. 

and 4.0 are too low; the results for these data points are not included in any 
of the fits. 

Table |^ summarises the coefficients extracted from various fits. A 'good 
fit' is defined in the same way as for the propagator fits: Q > 0.1. The fitting 
range for which Q falls below this cut-off is taken to signify when higher 
terms can be resolved. Thus, we do not present fits including an additional 
term over a smaller fitting range. The decay constant is only consistent 
with a linear dependence on 1/Mps for the first three data points, for which 
l/Mps^O.3. The quadratic term is needed for l/Mps'^0.6, i.e. there are 
significant 0(1/M^) contributions to the decay constant in the region of 
the B meson {l/{aMps) ~ 0.4). Since only the 0{1/M) contributions are 
correctly included in this analysis it is necessary to go to the next order in 
the NRQCD expansion to calculate fp accurately. A fit to all the data points 
is possible including the cubic term. 

A quadratic fit to the heaviest six data points is shown in figure ^ The 
extrapolation of the NRQCD results at finite heavy quark mass is clearly 
consistent with the static data point. Since the static theory is the heavy 
quark mass limit of NRQCD this is expected and we assume this is not spoilt 
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Order 


fit range 


n.d.o.f. 


Q 


a3/2Co 








1 


1-3 


1 


0.3 


0.324(6) 


-0.371(22) 






1-4 


2 


0.0 


0.317(6) 


-0.329(15) 






2 


1-4 


1 


0.7 


0.342(10) 


-0.561(71) 


0.458(137) 




1-5 


2 


0.7 


0.338(10) 


-0.525(55) 


0.379(75) 




1-6 


3 


0.7 


0.335(9) 


-0.486(40) 


0.304(60) 




1-7 


4 


0.7 


0.332(7) 


-0.456(23) 


0.254(26) 




3 


1-8 


4 


0.8 


0.346(9) 


-0.623(52) 


0.740(119) 


-0.425(104) 


1-9 


5 


0.9 


0.346(8) 


-0.614(45) 


0.716(100) 


-0.406(69) 



Table 7: The coefficients extracted from fits to a^/^(/v^M)p5 as a function 
of XjiaMps) for aMo = 1.0 and = 0.1585. Note tliat tlie renormalisation 
factors are not included. The fit range 1—3 denotes a fit to the data at 
aMo = 10.0, 7.0 and 3.0. The data at oMq = 4.0 and 3.5 are not included in 
the fit. 



by the inclusion of the renormalisation factors. A striking feature of the plot 
is the very large slope and the corresponding large deviations from the static 
limit of ~ 50% around l/M^. 

Note that with only the last two data points requiring a cubic fit, the 
coefficient C\ is not as well determined as for the quadratic fit. Hence, we 
only consider the quadratic fits and find a^/^Co = 0.34(1). Similarly, o?l'^C\ ~ 
—0.46(5). Thus, we arrive at acp = aCi/Co ~ —1.35(15), much larger than 
the naive expectation of 0{aAQCD) ~ —0.2. Converting to physical units, 
cp ~ —2.8(5) GeV, using = 1.8—2.4. We found the slope parameters do 
not depend significantly on the light quark mass. 

This slope is significantly larger than that found in previous lattice simu- 
lations which found cp 0.8 GeV |g and -1.14 GeV [0]. We believe the 



discrepancy is due to previous results being obtained by extrapolation from 
around the D meson with too small a range of masses when using Wilson or 
Clover heavy quarks; consistency between the static limit and results in this 
region has not been shown. In addition, our results indicate the deviations 
from the static limit may be so large at Mp, that a heavy quark expansion is 
not valid in this region. 

Onogi et al, using the Fermilab approach to an order equivalent to 0(1/M) 
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in NRQCD, also find a slope around —1 GeV data points around the 
B and D meson are used. The pseudoscalar renormalisation factors for this 
method have not yet been calculated and a proper comparison of the extrap- 
olation of the results with the static limit cannot yet be made. Consistency 
should be found between this method and NRQCD, and it is important to 
confirm the slope in the large mass region, Mps ^ Mb- 

To examine the large slope of the decay constant further we obtained the 
individual contributions to the pseudoscalar and vector matrix elements from 
each of the 0(1/M) interactions. Figure | shows (fy/M), {f\/M)pg and 
{f\/M)^g°^^ as a function of 1/Mps- From section |^ the slope of the spin- 
average decay constant depends only on the kinetic energy of the heavy quark, 
while the other two quantities include contributions from the spin dependent 
interactions. It is clear from the large mass region where the slope is linear 
that Gkin is much greater than Ghyp or Gcorr, and is the source of the large 
slope. Quahtatively, the spin dependent terms are expected to be suppressed 
since they break spin as well as flavour symmetry. Performing correlated 

fits to (fy/M), detailed in table ||, we find aGfci„ = Ci/Cq ~ —1.26(15) or 
-2.6(5) GeV. 

This result agrees with Hashimoto 0, who calculates the pseudoscalar 
decay constant at /3 = 6.0 in the quenched approximation and finds the co- 
efficient of the slope at 0{1/M) roughly —2.4(1.1) GeV. Since the correction 
to the current is omitted in his calculation, this coefficient corresponds to 

Note that even though only 0(1/M) terms are included in the NRQCD 
action we see quadratic and cubic behaviour in the decay constant in figure ^ 
this is due to higher powers of the existing 0(1/M) terms. These higher order 
powers are dominated by powers of the kinetic energy terms as is clear from 
figure H, which shows that the spin-dependent and spin-independent slopes 
of fvM show little difference even at relatively small values of M. Also, 
the missing 0(1/M^) terms in the NRQCD action are all spin-dependent so 
they cannot contribute to the spin-independent slope at 0(1/M^). We would 
therefore expect that our spin-independent slope was correct even through 
0(1/M2) except that there are spin-independent matrix element corrections 
at this order. We do not know the size of these corrections, but it is unlikely 
that they would be larger than the spin-independent 1/M^ terms that are 
present in our calculation already. Later we will allow for this by taking a 
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0.1 ^ ^ ^ I ^ ^ ^ ^ I ^ ^ ^ ^ I ^ ^ ^ — 

0.0 0.3 0.4 0.6 0.1 

l/(aMps) 

Figure 8: The pseudoscalar decay constant both with (crosses) and with- 
out (pluses) the current correction and the spin-average decay constant (cir- 
cles) as a function of l/{aMps) for ni = 0.1585. The solid lines indicate 
cubic fits to all the data points. Note that the renormalisation factors are 
not included. 



rather pessimistic 100% error in our 1/M^ coefficient. 

To isolate the effects of Ghyp and Gcorr, consider the ratio of the pseu- 
doscalar and vector decay constants, with and without the current correc- 
tions, shown in figure ^. Both data sets are consistent with 1 in the static 
limit, as expected. Note that without the current corrections the ratio, which 
is determined solely by the hyperfine interaction, is greater than 1 indicating 
Ghyp is positive. We find acp ~ -1-0.11(1) or .23(4) GeV in physical units, 
and hence aGhyp = d'p/% ~ 0.014(1) or 0.029(5) GeV. Including S{fVM), 
the slope is a physical quantity and acp = —0.40(4) or —0.85(15) GeV. 

By subtracting Cp and Cp we estimate aG 

corr ~ — 0.77(8). However, the 
correction to the current can be extracted directly and more accurately as the 
intercept of 2MQ5{f\fM)ps/{fVMY^^'''' shown in figure |10|. From table I, 
we obtain consistent results, aGcorr = —0.675(7) or —1.42(20) GeV. As in 
the finite heavy quark mass case, Gcorr agrees with — -E^^ = —0.528(5) to 
within approximately 0{as). 

A comparison can be made with the predictions of QCD sum rules. The 
results obtained by Neubert using HQET, shown in table 0, are in good 
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a3/2(/VM) 


0.335(7) 


-0.43(5) 


/ (fV M)ps \uncorr 
^ (/VM)i' ^ 


0.996(5) 


0.11(1) 


/ ifVM)ps \tot 
^ (fVM)v ' 


1.009(5) 


-0.40(4) 




aCo 






-0.675(7) 


0.19(1) 



Table 8: The coefficients extracted from fits to combinations of the pseu- 
doscalar and vector decay constants as a function of l/{aMps) for ki — 
0.1585. The errors include the variation in the coefficients using different 
orders in the fit function. 



1.2 




0.0 0.2 0.4 0.6 0.1 

l/(aMps) 

Figure 9: The ratio of the pseudoscalar to vector decay constant both 
with (crosses) and without (pluses) the current correction as a function 
of l/{aMps) for ki = 0.1585. The solid lines indicate a linear fit to 
((/v^)p5/(/v^)y)"'' and a quadratic fit to {{f^/M)ps/if^/M)vy°\ 
where all data points are included in the fits. The error in the extrapolation 
to the static limit is indicated for both fits. Note that the renormahsation 
factors are not included. 
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0.50 

0.0 0.3 0.4 0.6 0.1 

l/(aMps) 



Figure 10: 2MQa5{fy/M)ps/U^Tps'"'^ vs l/{aMps) for ki = 0.1585. The 
solid line indicates a linear fit to the six heaviest data points. Note that the 
renormalisation factors are not included. 

agreement with our analysis. As noted previously, the unphysical parameters 
Ghyp and Gcorr differ; from table ^ the former differs in sign as well as in 
magnitude with our result. It remains unclear, however, why the slope, 
and thus G^m; is much larger than naively expected; a similar analysis of the 
spectrum [0] showed deviations from the static limit in agreement with naive 
expectations. The decay constant may be unique in this respect, however, 
it is clear this should not be assumed, particularly for quantities such as 
semi-leptonic form factors of B mesons which are not protected by Luke's 
theorem. 

By chirally extrapolating the results and interpolating to l/{aMB), we 
find fs ~ 126-166 MeV. The static renormalisation factor, Z'^^^ = 0.73 Q, 
is used, where we guess aq*z^ = 2.0 at aMo = oo. To estimate the error in 
/b due to the truncation of the NRQCD series we assume a 100% error 
in the quadratic coefficients of the heavy quark mass expansion. Thus, we 
find a ~ 20% uncertainty in Jb, which is comparable to the uncertainties 
arising from the light quark sector. This is likely to be an overestimate; 

^This renormalisation constant corresponds to ZZcont in the notation used in this 
reference. Also, the anomalous dimension factor is set to 1. 
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Sum Rules 


NRQCD 


Cp 


-2.9(5) 


-2.8(5) 


III 

Cp 


-0.9(1) 


-0.85(15) 




-2.3(4) 


-2.6(5) 


Ghyp 


-0.07(3) 


0.029(5) 


Grcorr 


-0.50(7) 


-1.42(20) 



Table 9: A comparison in GeV of the decay constant slope parameters and 



coefficients extracted in this simulation with those computed by Neubert [19 



using HQET and QCD sum rules. The errors are dominated by the uncer- 
tainty in a^^. 



since the contributions from the kinetic energy term in the NRQCD action 
dominate the slope, the corrections to the matrix element at 0(1/M^) are 
unlikely to be larger than the quadratic contributions already present. Sim- 
ilarly, from the size of Gkin the next order spin-independent terms in the 
action, —p'^/8M^, are expected to lead to larger contributions to fs than 
naively predicted, but this is unlikely to be larger than the lower order con- 
tributions. The 0{as) ~ 20% error in the slope introduced by using the 
static renormalisation factor is less significant. In the static limit, we obtain 
fptat ^ 237-295 MeV at ki = 0.1585 which is close to k,. 

Our determination of Jb is consistent with the current world average 
of 200(40) MeV [^. Our analysis is being performed in conjunction with 
a simulation using (3 = 6.0 quenched configurations With such large 



systematic errors it is not possible to discern any difference between the 
partial inclusion of dynamical quarks in this analysis and quenched results. 
However, the decay constant is related to the wavefunction at the origin. 



JbVMb = vT2|^(0)|, and thus it is expected to be sensitive to short dis- 
tance physics and the quenched approximation. Booth |2^ has estimated 



the effect of quenching to be 10 — 15% in Jb and ~ 5% in Jbs/ fs using 
chiral perturbation theory coupled with HQET. However, the prediction is 
fairly dependent on the light quark mass. In the future we aim to reduce the 
systematic errors in both our quenched and dynamical simulations by using 
tadpole improved clover light fermions and including higher orders terms in 
the NRQCD expansion. Thus, we expect to isolate and remove the full effects 
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of quenching. 

With an experimental measurement of mixing also in the B'^—Bg sys- 
tem only the ratio fBsBsJlBBB is required to extract |Vts/Vtd|. Some 
cancellation of systematic errors is expected for this ratio, and hence it 
should be determined more accurately. Table [l^ presents the results for 
((/-\/M)s/(/A/M)d)*°*, where s and d denote results for a light quark mass 
extrapolated to the strange quark mass and the chiral limit respectively. The 
expected 0{ms/M) dependence on the heavy quark mass is not resolved in 
our results and the ratio is independent of Mq. The error indicated is purely 
statistical. There is an additional ~ la error due to the uncertainty in Kg. 



aMo 


0.8 


1.0 


1.2 


1.7 


2.0 


2.5 


({fy/M)s\tot 


1.26(3) 


1.26(3) 


1.25(3) 


1.25(3) 


1.25(3) 


1.26(4) 


aMo 


3.0 


3.5 


4.0 


7.0 


10.0 


oo 


(ifVM)s\tot 


1.26(4) 


1.26(4) 


1.25(4) 


1.24(5) 


1.25(3) 





Table 10: The ratio of the pseudoscalar decay constant for ki = Kg and 
Ki = He for all aMo, where Kg = 0.1577 is obtained using the ratio M^/Mp. 
The errors shown are purely statistical. 



6 Conclusions 

We presented a lattice study, partially including the effects of dynamical 
quarks, of the heavy-hght pseudoscalar and vector decay constants using 
Wilson light quarks and NRQCD heavy quarks. All 0{1/M) terms in the 
NRQCD action and matrix elements were included; an 0{aAQCD/M) wave- 
function renormalization induced by the evolution equation is the same order 
as the light quark discretization errors. We find consistency between the 
extrapolation of the NRQCD result and the static expected, although 

the static result is not very well determined. There are significant 0(1/M^) 
corrections to fVM around Mb, and this leads to a large systematic error in 
extracting /b; the use of Wilson light quarks introduces a systematic error 
of roughly the same magnitude. It is necessary to go to the next order in 
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NRQCD and use an improved light quark action in order to reliably extract 
fs. 

We found the slope of the decay constants with 1/M to be significantly 
larger than that found in previous lattice simulations. We believe this is 
because we are closer to the static limit. In particular, for the first time 
wc obtained the three 0{1/M) contributions to the linear component of the 
slope separately. We found that Gkin dominates and thus is responsible 
for the large slope. Good agreement is found between our results for Gkin 
and physical combinations of Ghyp and Gcorr and a prediction of QCD sum 
rules. The results show that the kinetic energy of the heavy quark gives rise 
to contributions to j\fM much greater than 0(Aqc'£)), contrary to naive 
expectations. 
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